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Introduction 

Approximately  2.5  million  Americans  live  with  epilepsy  and  epilepsy-related  deficits  today,  more 
than  disabled  by  Parkinson  disease  or  brain  tumors.  The  impact  of  epilepsy  in  the  US  is 
significant  with  a  total  cost  to  the  nation  for  seizures  and  epilepsy  of  approximately  $12.5  billion. 
Epilepsy  consists  of  more  than  40  clinical  syndromes  affecting  40  million  people  worldwide. 
Approximately  25  percent  of  individuals  receiving  antiepileptic  medication  have  inadequate 
seizure  control;  however,  80%  individuals  with  medication  resistant  epilepsy  might  be  cured 
through  surgery  if  one  were  able  to  precisely  localize  the  seizure  focus.  The  proposed  research 
will  significantly  advance  our  ability  to  localize  such  foci,  and  thereby  offer  curative  epilepsy 
surgery  for  this  devastating  disease.  Photoacoustic  tomography  (PAT)  uniquely  combines  the 
high  contrast  advantage  of  optical  imaging  and  the  high  resolution  advantage  of  ultrasound 
imaging  in  a  single  modality.  In  addition  to  high  resolution  structural  information,  the  proposed 
PAT  is  also  able  to  provide  functional  information  that  are  strongly  correlated  with  regional  or 
focal  seizure  activity,  including  blood  volume  and  blood  oxygenation  because  of  the  high 
sensitivity  of  optical  contrast  to  oxyhemoglobin  and  deoxyhemoglobin  concentrations.  The 
hypothesis  of  the  proposed  research  is  that  PAT  offers  the  possibility  to  non-invasively  track 
dynamical  changes  during  seizure  occurrence.  The  overall  goal  of  this  research  is  to  advance  a 
finite  element  based  photoacoustic  tomography  method  for  epilepsy  imaging,  using  both 
laboratory  and  in  vivo  experiments.  Specifically,  in  this  project  we  propose:  (1)  To  design, 
construct  and  test  a  transducer  array  system  for  both  2D  and  3D  PAT  imaging;  (2)  To  advance 
reconstruction  algorithms  and  associated  image  enhancement  schemes  for  quantitative  PAT;  (3) 
To  evaluate  and  optimize  the  integrated  functioning  of  the  hardware  and  software  components 
of  the  transducer  array-based  system,  using  simulation  and  phantom  experiments;  (4)  To  test 
and  validate  the  PAT  system  using  a  well  established  animal  model  of  temporal  lobe  epilepsy. 

Body 

This  report  describes  work  accomplished  in  Year  1  (Months  0-12)  of  the  project.  As  outlined  in 
the  approved  Statement  of  Work  (SOW),  the  tasks  during  this  period  of  time  include:  Task  1 
(Months  0-24):  Purchase  and  calibrate  Ti:Sapphire  tunable  laser;  design  and  test  PVDF 
transducers;  design  and  test  data  acquisition  subsystem;  assemble  the  entire  photoacoustic 
tomography  (PAT)  system;  Task  3  (Months  0-24):  Implement  reconstruction  codes  and 
associated  enhancing  schemes  including  dynamic  dual  meshing,  total  variation,  advanced 
regularization  techniques,  adaptive  meshing,  initial  parameter  optimization  and  reconstruction  of 
absolute  optical  absorption  coefficient  for  quantitative  high  resolution  functional  PAT. 

The  sections  below  consist  of  (1)  hardware  development  and  (2)  software  development  that 
reflect  the  tasks  associated  with  the  SOW  during  Months  0-12. 


1.  Hardware  Development  (Task  1) 

We  have  completed  the  design  of  the  proposed  array  based  PAT  system  based  on  the 
extensive  modeling  and  testing  of  the  key  components  needed  for  the  proposed  system  using 
the  existing  single  transducer  scanning  system.  Purchases  of  all  these  components  were 
already  out,  and  we  should  be  in  a  position  to  test  each  of  the  received  components  in  next  2-3 
months  and  ready  to  assemble  the  entire  system  during  the  summer.  We  are  confident  that  we 
will  be  able  to  build  an  operational  PAT  system  by  the  end  of  this  year. 
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Here  we  detail  our  designing  work  for  an  array  based  PAT  system  which  is  schematically  shown 
in  Fig.  la. 


Fig.  1  (a)  Schematic  of  the  PAT  system,  (b)  and  (c):  3D  display  of  the 
system/measurement  chamber  with  (b)  and  without  (c)  the  top  covered. 

1.1  Excitation  light  source 

As  proposed,  a  tunable  Ti:Sapphire  laser  (690-900nm,  8-25ns  pulse  width,  120mJ  between  760 
and  800nm)  is  already  in  place  that  will  be  used  as  the  light  source. 

1.2  Transducer  array 

Originally  we  proposed  to  customize  a  PVDF-based  ring-type  phase  array  with  three  rows  of 
transducer  elements.  However,  our  recent  study  showed  that  the  sensitivity  of  PVDF-based 
array  is  not  good  enough  to  guarantee  high  quality  data  collection  especially  when  the  working 
wavelength  moves  to  near-infrared  (NIR)  (i.e.,  690-900nm).  So  we  decided  considering  to 
employ  a  phased  array  made  of  piezocomposite  materials  instead  of  PVDF  for  our  PAT  system. 
The  cost  of  piezocomposite  ring-type  phased  array  (>  $70,000),  is  much  higher  than  a  PVDF- 
based  array.  In  order  to  reduce  the  cost  for  a  phased  transducer  array,  we  have  conducted 
numerical  modeling  and  phantom  experiments  to  see  if  we  can  use  a  flat  1 .5D  phased  array  to 
obtain  imaging  quality  that  is  comparable  to  a  ring-type  phased  array. 

In  the  numerical  simulations,  a  finite  element  code  in  time  domain  was  used  to  generate  the 
ultrasound  signals  at  all  detector  locations  and  the  delay-and-sum  algorithm  was  used  to 
reconstruct  the  photoacoustic  images.  Here  the  images  obtained  from  a  ring-type  array  and  an 
array  formed  by  8-piece  flat  elements  are  compared.  Two  kinds  of  targets,  circles  and  crosses, 
were  embedded  in  a  background  medium.  Fig.  2  shows  the  exact  (a  &  d)  and  reconstructed 
photoacoustic  images  (b,  c,  e,  &  f)  when  each  element  was  treated  as  an  ideal  point  transducer 
that  can  receive  signals  coming  from  all  angles.  We  see  that  the  reconstructed  images  from  the 
two  types  of  arrays  are  almost  identical. 
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Fig.  2.  Simulation  results  without  applying  transducer  directivity,  (a)  &  (d):  the  exact 
images,  (b)  &  (e):  reconstructed  images  based  on  a  ring-type  array,  (c)  &  (f): 
reconstructed  images  an  8-piece  flat  element-based  array. 

In  a  real  situation,  each  transducer  element  has  directivity,  meaning  that  it  can  only  collect 
ultrasound  signal  from  a  limited  angle,  depending  on  its  size  and  central  frequency.  Fig.  3 
shows  the  typical  2D  fields  from  a  5MHz  transducer  with  2mm  in  diameter  and  a  4MHz 
transducer  with  1 .5mm  in  diameter.  Generally,  a  transducer  with  smaller  size  and  lower  central 
frequency  has  bigger  receiving  angle,  but  lower  sensitivity  and  poorer  spatial  resolution. 
Considering  all  the  trade-offs  among  spatial  resolution,  sensitivity  and  receiving  angle, 
transducer  elements  with  a  central  frequency  of  4MHz  and  a  dimension  of  1.5mmx2.0mm  will 
be  used  in  our  application.  Fig.  4  gives  the  field  coverage  of  the  ring-type  array  and  8-piece  flat 
array. 


Fig.  3.  Calculated  2D  fields  from  (a)  5MHz 
<t>2mm  and  (b)  4MHz  <D1 ,0mm  transducer 
element. 


Fig.  4.  2D  field  coverage  of  (a)  the  ring-type  array 
and  (b)  the  8-piece  flat  array.  The  dimension  of  the 
transducer  elements  is  1 .5x2. 0mm  and  the  central 
frequency  is  4MHz. 


When  the  directivity  of  each  transducer  element  is  considered  in  the  numerical  simulations,  the 
results  show  that  the  8-piece  flat  array  can  still  provide  satisfactory  image  quality  (Fig.  5). 
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Fig.  5.  Simulation  results  with  directivity  considered,  (a)  &  (d):  the  exact  images,  (b)  & 
(e):  reconstructed  images  based  on  a  single  ring-type  array,  (c)  &  (f)  reconstructed 
images  based  on  8-piece  flat  array. 


We  have  conducted  phantom  experiments  to  confirm  the  simulation  findings.  In  the  experiments, 
multiple  hairs  were  embedded  in  a  gel  phantom  background  (see  a  photograph  shown  in  Fig. 
6a).  A  single  transducer  with  2mm  in  diameter,  driven  by  a  motorized  rotator  and  a  linear  stage 
was  used  to  collect  the  ultrasound  signals  along  a  circular  path  (simulating  a  ring-type  array) 
and  along  an  arc-straight  line  path  (simulating  an  8-,  16-,  and  32-piece  flat  array).  The 
experimental  images  are  shown  in  Fig.  6  (b-e).  Taking  a  close  look  at  the  recovered  images,  we 
found  that  some  structures  were  missed  in  the  image  from  a  simulated  8-piece  flat  array.  We 
realized  that  this  was  due  to  the  missing  of  ultrasound  signals  from  some  directions  for  this 
array  configuration.  We  note  that  the  image  from  a  simulated  16-piece  flat  array  gives  almost 
the  same  quality  as  the  simulated  ring-type  array.  Based  on  these  phantom  experiments,  we 
decided  that  al 6-piece  flat  array  (4x3  elements  each  piece)  will  be  used  for  our  PAT  system. 
The  cost  for  this  array  is  reduced  by  >  $12,000  compared  to  a  ring-type  array. 


-ID  0  ID  -ID  0  ID  -ID  0  ID 

X  (mmj  X  (mm)  X  (mm) 


Fig.  6.  Experimental  results,  (a)  A  photograph  of  the  hair-gel  phantom.  The 
recovered  imaged  from  the  simulated  ring-type  array  (b),  8-piece  flat  array 
(c),  16-piece  array  (d)  and  32-piece  array  (e). 


4 


Huabei  Jiang,  PhD/Annual  Report 


The  16-piece  flat  transducer  arrays  will  be  arranged  along  a  circular  path  with  a  radius  of  41mm. 
We  asked  quotations  for  such  a  flat  array  from  three  companies  including  GE  (Lewistown,  PA), 
Olympus-IMS  (Waltham,  MA),  and  Imasonic  SAS  (I'Ognon,  France),  and  we  decided  to  buy 
from  Imasonic  SAS  as  it  is  the  most  cost-effective.  In  this  array,  each  piece  comes  with  3  (row) 
x  4  (column)  elements,  which  are  directly  cabled,  and  each  of  them  acts  as  an  independent 
detecting  channel,  making  a  total  of  192  channels  on  three  rows  for  the  whole  array  subsystem. 
The  central  frequency  of  the  elements  is  4MHz  with  60%  bandwidth.  The  dimension  of  each 
element  is  2.0  x  1.5mm.  Fig.  7  shows  the  structure  of  the  detection  unit  as  well  as  the 
distribution  of  the  elements  in  each  transducer  array. 


Fig.  7.  The  arrangement  of  16-piece  flat  array.  The  insert  shows  the  distribution  of  elements 
in  one  array. 


Fig.  8.  PCB  under  testing  for  Fig.  9.  Testing  result  of  ultrasound  signal:  red, 

12-channel  preamplifiers.  from  the  Olympus  5058PR  ;  blue:  from  the 

PCIAD850  1-8  channel,  respectively. 


1.3  Data  acquisition  (Preamplifiers,  multiplexer  and  A/D  board) 

The  originally  proposed  data  acquisition  system,  DiPhAS  by  Fraunhofer  Institute  for  Biomedical 
Engineering  (St.  Ingbert,  Germany),  contains  preamplifiers,  signal  processing  and  multi-channel 
A/D  boards  in  one  package  and  provides  quality  needed  for  the  proposed  study.  However,  the 
price  of  this  system  is  increased  significantly  since  early  2009  (before  the  start  of  this  project), 
and  the  approved  budget  of  the  project  does  not  allow  us  to  buy  this  system.  Since  there  is  no 
commercial  one-package  data  acquisition  affordable  to  us  in  the  current  market,  we  decided  to 
do  the  following  to  resolve  the  issue:  purchase  the  A/D  boards  and  build  preamplifiers  ourselves. 


5 


Huabei  Jiang,  PhD/Annual  Report 


We  have  a  total  of  192  signal  channels  (corresponding  to  192  transducer  elements  from  the 
array).  Each  channel  will  need  a  preamplifier.  The  design  of  preamplifier  is  based  on  ultralow 
distortion,  ultralow  noise  and  high  speed  AD8099  operational  amplifier  from  Analog  Device, 
which  can  provide  26dB  gain.  Bandwidth  filters  have  been  added  to  the  preamplifier  circuits  to 
further  decrease  the  noise  level.  Fig.  8  is  the  designed  12-channel  PCB  (printed  circuit  board).  A 
total  of  sixteen  12-channel  PCBs  will  be  used  and  placed  right  beside  the  16-piece  flat 
transducer  array,  see  Fig.  1c. 

Ideally  we  would  need  a  192-channel  A/D  boards,  but  the  cost  is  far  beyond  the  approved 
budget  for  equipment.  Instead  we  just  can  afford  a  64-channel  A/D  boards  (US  Ultratek  Concord, 
CA).  The  data  acquisition  from  all  192  channels  is  realized  via  an  integration  of  the  64-channel 
boards  and  a  multiplexer.  The  64-channel  A/D  boards,  PCIAD850/PCIAD1650,  can  provide 
50MHz  sampling  rate  for  all  channels  at  a  resolution  up  to  1 0bit.  The  on-board  gain  can  be 
adjusted  separated  for  each  channel  from  -12dB  to  84dB.  A  demo  board  from  this  manufacturer 
has  been  tested  and  has  been  compared  to  the  5058PR,  a  commercial  single  channel 
pulser/receiver  from  Olympus.  The  comparative  measurements  are  shown  in  Fig.  9,  indicating 
that  the  performance  of  these  data  acquisition  boards  is  satisfactory  for  our  project. 

As  indicated  above,  the  ultrasound  signals  from  the  192  channels  will  be  selected  and  delivered 
to  64  parallel  data  acquisition  channels  through  a  multiplexer  placed  under  the  imaging  interface 
as  shown  in  Fig.  Ib-c.  The  multiplexer  will  be  purchased  from  Cytec  Corp.  (Penfield,  NY).  We 
are  also  making  effort  on  another  option  by  integrating  multiplexer  into  the  preamplifier  PCBs. 
The  advantages  of  this  solution  include  the  more  compactness  of  the  system,  lower  noise  and 
lower  cost. 

The  computer,  which  is  used  for  system  control  and  data  storage,  will  be  an  industrial  4U  PC 
with  a  dual-core  processor  and  at  least  8  PCI  slots.  We  have  a  solution  from  Acnodes 
Corporation  (Walnut,  CA).  A  LabVIEW  program  will  be  used  to  control  the  whole  system. 

2.  Software  Development  (Task  3) 

In  Task  3,  we  will  implement  6  image  enhancement  schemes.  In  Year  1,  we  have  chosen  to 
complete  the  implementation  of  3  schemes  including  total  variation  (TV),  advanced 
regularization  techniques,  and  adaptive  meshing.  We  will  implement  the  remaining  3  schemes 
in  Year  2. 

2.1  Total  variation  minimization 

The  concept  of  TV  was  originally  conceived  as  a  way  of  restoring/enhancing  images  (e.g., 
denoising,  deblurring  etc.).  Conventional  PAT  reconstruction  algorithms  are  based  on  a 
nonlinear  least  squares  criterion  which  stands  on  the  statistical  argument  that  the  least  squares 
estimation  is  the  best  over  an  entire  ensemble  of  all  possible  pictures.  TV,  on  the  other  hand, 
measures  the  oscillations  of  a  given  objective  function  and  does  not  unduly  punish 
discontinuities.  Hence,  we  have  hypothesized  that  a  hybrid  of  these  two  minimization  schemes 
should  be  able  to  provide  higher  quality  image  reconstruction.  We  have  successfully 
implemented  the  TV  scheme  into  the  existing  reconstruction  codes  (see  the  Appendix  for 
mathematical  detail).  We  have  evaluated  the  TV  scheme  using  experimental  data.  For 
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comparative  purposes,  reconstructions  without  the  total-variation-minimization  enhancement  are 
also  presented. 

Three  phantom  experiments  were  conducted.  In  the  first  two  experiments,  we  embedded  one  or 
two  objects  with  a  size  ranging  from  3  to  0.5  mm  in  a  50  mm-diameter  solid  cylindrical  phantom. 
The  phantom  materials  used  consisted  of  Intralipid  as  scatterer  and  India  ink  as  absorber  with 
Agar  powder  (1-2%)  for  solidifying  the  Intralipid  and  India  ink  solution.  The  absorption  coefficient 
of  the  background  phantom  was  0.01  mm'1,  while  the  absorption  coefficient  of  the  target(s)  was 
0.03  mm'1.  In  the  last  experiment,  we  used  a  single-target-containing  phantom,  aiming  to  test 
the  capability  of  detecting  target  having  low  optical  contrasts  relative  to  the  background 
phantom.  In  this  case,  the  target  had  an  absorption  coefficient  of  0.015  mm'1.  The  reduced 
scattering  coefficients  of  the  background  phantom  and  targets  were  1.0  and  3.0  mm'1  for  the 
first  two  experiments,  and  1.0  and  2.0  mm'1  for  the  last  experiment.  In  these  experiments,  a  total 
of  120  receivers  were  equally  distributed  along  the  surface  of  the  circular  background  region. 

Figure  10  shows  the  reconstructed  photoacoustic  images  for  all  the  three  experimental  cases, 
while  Figure  1 1  presents  quantitative  absorption  profiles  along  transacts  through  one  target  for 
the  images  shown  in  Fig.  10.  We  see  that  considerably  enhanced  images  are  achieved  with  the 
total-variation  minimization,  especially  when  the  target  is  tiny  (case  2),  or  the  contrast  level 
between  the  target  and  the  background  is  low  (case  3). 


Fig.  10  Reconstructed  absorbed  energy  density  images,  (a)  case  1  without  TV  minimization,  (b)  case  2  without 
TV  minimization,  (c)  case  3  without  TV  minimization,  (d)  case  1  with  TV  minimization,  (e)  case  2  with  TV 
minimization,  (f)  case  3  with  TV  minimization. 
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X  position  (mm) 


Fig.  11  Recovered  absorbed  energy  density  profiles  along  (a)  y=-7.0mm  for  case  1  (3mm  diameter  target),  (b) 
y=8.0mm  for  case  1  (2mm  diameter  target),  (c)  y=1 ,0mm  for  case  2,  (d)  y=6.5mm  for  case  3. 

2.2  Advanced  regularization  techniques 

Image  reconstruction  is  typically  an  ill-posed  problem.  That  is,  its  solution  procedures  can 
become  singular  or  nearly-singular  such  that  small  errors  in  the  measured  data  are  magnified 
significantly  in  the  resulting  image.  Regularization  techniques  are  practical  methods  which  have 
proven  to  be  able  to  combat  the  ill-posed  problems  associated  with  noisy  data.  We  have 
implemented  a  prototype  hybrid  technique  which  synthesizes  standard  Marquardt  and  Tikhonov 
regularization  schemes  in  the  existing  reconstruction  codes.  As  proposed,  we  have  attempted  to 
further  evaluate  and  optimize  this  method  in  Year  1.  Fig.  12  shows  a  set  of  typical 
reconstructions  for  the  optimization  study  where  the  regularization  parameter,  A,  is  varied  from 
0.6><102to  0.6x10'6  for  a  test  case  having  three  targets  embedded  in  a  circular  background 
medium.  We  note  that  these  images  are  similar  qualitatively;  however,  quantitatively  in  terms  of 
the  recovered  optical  absorption  (see  the  color  bar  on  the  right),  A=0.6x10_1(d)  and  A=0.6x10‘2 
(e)  give  the  best  results. 

We  have  also  implemented  the  weighted  least  squares  criterion  in  our  reconstruction  algorithms. 
Fig.  13  gives  the  reconstructed  images  from  a  typical  test.  We  see  that  both  images  are  almost 
identical,  but  as  expected,  the  reconstruction  with  the  weighted  least  squares  criterion 
converged  50%  faster  than  that  without  such  criterion. 

We  have  implemented  a  spatial  low  pass  filter  in  our  algorithms,  which  acts  on  the 
photoacoustic  property  distribution  during  the  iterative  reconstruction  process.  Fig.  14  shows  the 
images  from  the  three-target  example.  We  note  that  this  filtering  is  useful  in  reducing  the  noise 
effect  in  the  background  and  improved  the  quantitative  nature  of  the  recovered  photoacoustic 
property  (see  the  color  bar  on  the  right). 
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Fig.  12.  Reconstructed  photoacoustic  images.  Left  target:  5mm  in  diameter  and  4:1  contrast;  middle  target:  3  mm 
in  diameter  and  3:1  contrast;  right  target:  8  mm  in  diameter  and  5:1  contrast,  (a):  A=0.6x102.  (b):  A=0.6x101.  (c): 
A=0.6.  (d):  A=0.6x10_1.  (e):  A=0.6x10-2.  (f):  A  =0.6x1  O'4,  (g):  A  =0.6x1  O'6. 


Fig.  13.  Reconstructed  photoacoustic  images  without  (left)  and  with  (right)  the  weighted  least  squares  criterion. 


Fig.  14.  Reconstructed  images  without  (left)  and  with  (right)  spatial  filtering. 

2.3  Adaptive  meshing 

We  have  successfully  implemented  the  adaptive  meshing  scheme  in  the  reconstruction  codes. 
In  this  scheme,  we  locally  refine  the  property  mesh  in  areas  where  the  interelement  property 
gradient  exceeds  a  preset  maximum  value.  At  boundaries  of  heterogeneities  one  would  expect 


9 


Huabei  Jiang,  PhD/Annual  Report 


large  gradients  in  the  property  distribution;  hence  this  scheme  will  adaptively  refine  such  areas. 
Reconstructed  results  from  two  representative  cases  are  shown  in  Fig.  15  where  the  uniform 
and  adaptive  meshes  used  are  also  given.  We  see  that  the  target  is  clearly  better  recovered 
with  the  adaptive  mesh  (d  &  f)  in  terms  of  its  shape  and  size  over  that  with  the  uniform  mesh  (c 
&  e)  for  both  cases.  We  also  note  that  the  background  region  for  both  cases  is  more  smoothly 
recovered  with  the  adaptive  mesh  compared  to  that  with  the  uniform  mesh. 
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Fig.  15.  Reconstructed  photoacoustic  images  from  experimental  data,  (a):  the  uniform  finite  element  mesh  used 
(2458  nodes),  (b):  the  mesh  locally  adapted  around  the  target  region  (3882  nodes),  (c):  the  image  for  case  1  with 
uniform  mesh,  (d):  the  image  for  case  1  with  adaptive  mesh,  (e):  the  image  for  case  2  with  uniform  mesh,  (d):  the 
image  for  case  2  with  adaptive  mesh.  Case  1:  a  2-mm-diameter  target  having  1.5:1  contrast.  Case  2:  a  2-mm- 
diameter  target  having  2:1  contrast. 


Key  Research  Accomplishments 

1.  We  have  completed  the  design  for  an  array  based  PAT  system,  which  is  the  most  important 
part  of  Task  1. 

2.  We  have  developed  three  novel  schemes  that  can  enhance  our  current  reconstruction 
software. 

3.  We  have  conducted  phantom  experiments  that  confirmed  our  software  enhancement. 


Reportable  Outcomes  (see  the  Appendix  to  this  Summary  Report) 
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Conclusions 

We  have  made  a  significant  progress  that  has  fulfilled  the  statement  of  work  proposed  for  Year 
1  of  this  project.  Given  the  successful  first  year,  we  will  be  able  to  fulfill  or  exceed  the  work 
statement  for  Year  2. 

Appendix 

L.  Yao  and  H.  Jiang,  Enhanced  time-domain  photoacoustic  tomography  through  total-variation 
minimization,  in  Biomedical  Optics  2010  Technical  Digest  (Optical  Society  of  America,  Miami, 
Florida,  2010). 
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Abstract:  A  total-variation-minimization-based  iterative  algorithm  is  describ  ed  in  this  paper  that 
enhances  the  quality  of  reconstructed  images  with  time-domain  data  over  that  obtained  previously 
with  a  regularized  least-squares  approach. 
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1.  Introduction 

Photoacoustic  tomography  (PAT)  is  an  emerging  non-invasive  imaging  technique  that  combines  the  merits  of  high 
optical  contrast  and  high  ultrasound  resolution  in  a  single  modality  [1-3].  In  PAT,  the  Helmholtz-like  photoacoustic 
(PA)  wave  e  quation  has  been  c  ommonly  used  as  a  n  ace  urate  m  odel  for  desc  ribing  1  aser- induced  a  coustic  wa  ve 
propagation  i  n  t  issue.  While  anal  ytical  reco  nstruction  m  ethods  h  ave  bee  n  u  sed  f  or  photoacoustic  im  age 
reconstructions,  finite  element  method  (FEM)  based  method  appears  to  be  particularly  powerful  in  this  regard  [4,  5]. 
In  our  previous  w  ork,  e  xtensive  si  mulations  a  nd  e  xperiments  wi  th  f  requency-domain  m  easurements  have  been 
conducted,  and  maps  of  the  absorbed  energy  density  in  multicentimeter  phantom  geometries  has  been  achieved  with 
a  Newton  iterativ  e  algo  rithm  u  sed  i  n  conjunction  with  fi  nite  elem  ent  m  ethods  [4-6].  Recen  tly  th  ere  has  been 
considerable  interest  in  time-domain  schemes,  since  this  algorithm  shows  its  advantages  when  recovering  absorbed 
energy  de  nsity  im  ages,  sue  h  as  less  bac  kground  a  rtifacts,  m  ore  accurately  recovered  object  size  and  contra  st 
between  the  object  and  background  [7].  Although  encouraging  successes  with  model-based  image  reconstruction 
have  been  reported,  noise  can  be  a  significant  limitation  for  these  techniques.  Our  own  experience  in  this  regard  has 
revealed  that  the  quantitative  quality  of  reconstruction  is  affected  considerably  when  data  with  large  noise  level  is 
used. 

In  this  study,  we  propose  the  incorporation  of  a  total-variation-minimization  scheme  into  our  time-domain  PAT 
reconstruction  algorithm  to  reduce  the  effects  of  noise  and  to  enhance  the  quality  of  the  reconstructed  images.  We 
demonstrate  the  new  scheme  using  several  tissue-like  phantom  experiments. 

2.  Total- Variation-Minimization  Scheme 

To  describe  o  ur  t  otal-variation-minimization  m  ethod,  we  fi  rst  b  riefly  i  ntroduce  t  he  t  ime-domain  ph  otoacoustic 
reconstruction  al  gorithm  wi  th  the  re  gularized  Newton  m  ethod  we  have  use  d  previously  [4-  7].  The  t  ime-domain 
photoacoustic  wave  equation  in  tissue  can  be  described  as  follows 

WM (1) 

v0  ot  Cp  ot 

where  p  is  t  he  press  ure  wave;  v0  is  th  e  sp  eed  of  acou  stic  wav  e  in  th  e  medium,;  /3  is  th  e  t  hermal  exp  ansion 

coefficient;  C  is  the  specific  heat;  ®  is  the  absorbed  energy  density;  j{t )  =  S(t  ~tQ)  is  assumed  in  our  study. 

To  form  an  image  from  a  presumably  uniform  initial  guess  of  the  absorbed  energy  density  distribution  we  need  a 
method  of  updating  ®  from  its  starting  value.  This  update  is  accomplished  through  the  least-squares  minimization 
of 

HP,f)  =  UP°  -  P°s)  ■  ® 

7=1 

where  p°j  and  pCj  are  observed  and  computed  acoustic  field  data  for  i  =  1,2,  •  •  •  M  boundary  locations.  Using  a 
regularized  Newton  method,  we  obtained  the  following  equation  for  updating  ®  : 

+  i i)aI  =  Zt(p0-pc),  (3) 
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where  p  =  , p2  , •  •  •  pM  )  ,p  =  ^  , p2  ? *  *  * PM  )  > is  the  update  vector  for  the  absorbed  optical  energy 

density,  3  is  the  Jacobian  matrix  formed  by  dp  /  d ®  at  the  boundary  measurement  sites;  X  is  the  regularization 
parameter  determined  by  combined  Marquardt  and  Tikhonov  regularization  schemes,  and  I  is  the  identity  matrix. 

We  now  incorporate  the  total  variation  of  ®  as  penalty  term  by  defining  a  new  functional  [8-11]: 

F(p,o)=F(p,o)+z(o),  (4) 

where  =  j" ^ry^|V<t>  +  S2  dxdy  is  the  penalty  term,  COm  and  8  are  typically  positive  parameters  that  need 

to  be  det  ermined  numerically.  The  m  inimization  of  Eq.  (4)  proceeds  in  standard  fashion  by  the  differentiation  of 

F  with  respect  to  each  nodal  parameter  that  constitutes  the  ®  distribution;  simultaneously  all  these  relations  are  set 
to  zero.  Then  we  find  that  the  minimizers  of  Eq.  (4)  are  iteratively  found  by  solving  the  matrix  system 

(3r3  +  i?  +  2/^  =  3r(/?°-jpe)-F,  (6) 

where  V  is  formed  by  dL/ 3®  and  R  is  formed  by  dV / 5®  .  After  Eq.  (6)  has  been  reached,  our  solution  procedure 

(standard  Gausssian  elimination)  and  regularization  parameter  selection  (Marquardt  method)  are  identical  to  those 
used  previously.  Hence,  it  b  ecomes  clear  that  the  only  additions  to  our  new  algorithm  result  from  the  assembly  of 
matrix  R  and  the  construction  of  column  vector  V. 

3.  Results 

In  t  his  sectio  n  th  e  enh  anced  recon  struction  algo  rithm  w  ith  to  tal-variation  m  inimization  is  ev  aluated  with  so  me 
experiment  data.  F or  comparative  purposes,  reconstructions  without  the  total-variation-minimization  enhancement 
(i.e.,  those  which  use  our  previous  regularized  Newton  method)  are  also  presented. 

The  experimental  setup  used  for  collecting  experimental  data  is  a  pulsed  ND:YAG  laser  based  single  transducer 
(1MHz)  scanning  system,  which  was  described  in  detail  elsewhere  [5].  Three  phantom  experiments  were  conducted. 
In  the  first  two  experiments,  we  embedded  one  or  two  objects  with  a  si  ze  ranging  from  3  t  o  0.5  mm  in  a  5  0  mm- 
diameter  solid  cylindrical  phantom.  The  phantom  materials  used  consisted  of  Intralipid  as  scatterer  and  India  ink  as 
absorber  with  Agar  powder  (1-2%)  for  solidifying  the  Intralipid  and  India  ink  solution.  The  absorption  coefficient  of 
the  background  phantom  was  0.01  mm"1,  while  the  absorption  coefficient  of  the  target(s)  was  0.03  mm'1.  In  the  last 
experiment,  we  used  a  single-target-containing  phantom,  aiming  to  test  the  capability  of  detecting  target  having  low 
optical  contrasts  relative  to  the  background  phantom.  In  this  case,  the  target  had  an  absorption  coefficient  of  0.015 
mm"1.  The  reduced  scattering  coefficients  of  the  background  phantom  and  targets  were  1.0  and  3.0  mm"1  for  the  first 
two  experiments,  and  1.0  and  2.0  mm"1  for  the  last  experiment. 

A  total  of  120  receivers  were  equally  distributed  along  the  surface  of  the  circular  background  region.  We  chose 

T  =  3.0xl04^  as  the  tim  e  range  a  nd  At  =  200ns  as  th  e  tim  e  in  terval  fo  r  th  e  im  age  reconstructions  si  nee 
longer  tim  e  ran  ge  and  sho  rter  tim  e  in  terval  do  not  app  ear  to  provide  better  recon  struction  results.  We  also 
implemented  a  dual  meshing  method  for  fast  yet  accurate  inverse  computation.  The  fine  mesh  used  for  the  forward 
calculations  consisted  of  5977  nodes  and  1 1  712  elements,  while  the  coarse  mesh  used  for  the  inverse  calculations 
had  1525  nodes  and  2928  elements.  All  the  images  obtained  from  regularized  Newton  method  are  the  results  of  three 
iterations,  while  those  obtained  from  the  total-variation  minimization  are  results  of  fifteen  iterations.  (We  found  that 
larger  number  of  iteration  changed  the  solutions  only  by  less  than  0.5%).  Parallel  code  was  used  to  perform  these 
calculations  and  each  iteration  cost  -40  minutes  on  Beowulf  clusters  with  4  CPUs.  For  the  cases  presented,  values  of 

.0, 8  —  0.001  for  the  first  and  third  case,  and  CO 0  =2.0,8  =  0.001  for  the  second  case  appear  to  provide 
excellent  results. 

Figure  1  shows  the  reconstructed  absorption  coefficient  images  for  all  the  three  experimental  cases,  while  Figure 
2  presents  quantitative  absorption  coefficient  profiles  along  transacts  through  one  target  for  the  images  shown  in  Fig. 
1.  We  see  that  considerably  enhanced  images  are  achieved  with  the  total-variation  minimization,  especially  when  the 
target  is  tiny  (case  2),  or  the  contrast  level  between  the  target  and  the  background  is  low  (case  3). 
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(d)  (e)  (f) 

Fig.l  Reconstructed  absorbed  energy  density  images,  (a)  reconstruction  of  case  1  with  no  TV  minimization,  (b)  reconstruction  of  case  2 
with  no  TV  minimization,  (c)  reconstruction  of  case  3  with  no  TV  minimization,  (d)  reconstruction  of  case  1  with  TV  minimization,  (e) 
reconstruction  of  case  2  with  TV  minimization,  (f)  reconstruction  of  case  3  with  TV  minimization. 


(a) 


(b) 


(c) 


(d) 


Fig.  2  Recovered  absorbed  energy  density  profiles  along  (a)  y=-7.0mm  for  case  1  (3mm  diameter  target),  (b)  y=8.0mm  for  case  1  (2mm  diameter 

target),  (c)  y=l  .0mm  for  case  2,  (d)  y=6.5mm  for  case  3. 


In  summary,  we  have  demonstrated  time-domain  photoacoustic  image  reconstructions  using  an  algorithm  that  is 
based  on  total-variation  minimization.  The  results  have  shown  that  enhancements  of  the  reconstruction  images  have 
been  achie  ved  whe  n  t  he  ne  w  algorithm  is  com  pared  with  the  re  gularized  least-squa  res  m  inimization.  T  his  ne  w 
method  may  have  largest  impact  on  image  recovery  for  cases  having  low  contrast  levels  between  the  target  and  the 
background,  or  cases  having  tiny  heterogeneities. 
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